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Abstract 

In this article we propose a new symbolic-numeric algorithm to find 
positive equilibria of a n-dimensional dynamical system. This algorithm 
implies a symbolic manipulation of ODE in order to give a local approxi¬ 
mation of differential equations with power-law dynamics (S-systems). A 
numerical calculus is then needed to converge towards an equilibrium, giv¬ 
ing at the same time a S-system approximating the initial system around 
this equilibrium. This algorithm is applied to a real biological example in 
14 dimensions which is a subsystem of a metabolic pathway in Arabidopsis 
Thaliana. 


1 Introduction 

The modelling and study of biological or biochemical systems has become an 
exciting challenge in applied mathematics. The complexity of real biological 
dynamical systems lies essentially in the non-linearities of the dynamics as well 
as in the huge dimension of systems, often leading to a numerical approach. 
However, as the understanding of cellular mechanisms grows, it has become 
obvious that the modelling step strongly needs symbolic tools in order to ma¬ 
nipulate more and more information and data, and to improve computational 
tools. Therefore a new area emerged, called ’’systems biology”. It involves dif¬ 
ferent fields of applied mathematics, from computer algebra (see for instance 
[10]) to numerical computation ([7]). 

In the past decades, a lot of different frameworks have been developped to 
study behaviors of complex biochemical processes. Let us cite here three of 
them: the discrete networks (see the work of R. Thomas [15]), the piecewise lin¬ 
ear systems (the so-called Glass networks [6], see also [5]) and sigmoidal switch 
systems ([12]). The main goal of all these approaches is to propose a (more 
or less) generic class of dynamical systems, either discrete or differential, that 
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model some behaviors of complex interaction systems. Once this class is clearly 
defined, its mathematical relevance generally allows both theoretical and nu¬ 
merical analysis. 

The class of systems we use in this paper is the set of S-systems (see [1], [16], 
[17]). The basic idea of this model is to represent interactions between bio¬ 
chemical species with power-law dynamics. Their mathematical expression is 
quite general, but sufficiently simple to allow theoretical and practical investi¬ 
gations. We propose in this article a symbolic-numeric algorithm that is based 
upon S-systems theory. Its goal is to compute the positive equilibria of a n- 
dimensional system of ordinary differential equations (ODE). As it converges 
towards an equilibrium, it provides a S-system that approaches the original 
dynamics around this equilibrium. As we will see, the local approximation of 
some dynamics with power-laws can be made symbolically in any point of the 
phase space. It can also include treatment of symbolic parameters. However, 
iterating this process in order to converge towards equilibria needs a numerical 
computation, which prevents the use of pure symbolic tools to the end. 

In the following, we give a definition of the S-system class as it can be found 
in the litterature(see for instance [16]). We then propose a symbolic-numeric 
algorithm that computes an iteration leading to the positive equilibria of a dy¬ 
namical system. We will see an application of this algorithm on a biological 
example in dimension 14. We finally conclude with some remarks on our algo¬ 
rithm and some future works. 


2 S-systems 

2.1 Definition 

We give here a definition of the class of S-systems. 

Definition 2.1 A n-dimensional S-system S(a, (3, G, H) is a dynamical system 
defined by the n differential equations: 

n n 

Xi = Q-i xf* - (3i xf 13 , i = 1... n 

j =i i=i 

with a = (aq, ... ,a„) G (R+) n , f3 = (/3i ... (3 n ) G (R+) n 
and G = (9ij) i>j=1 _ n G M„(R), H = (h ;j ) ; j _ l n G -M„(R). 

R!j_ denotes the set of strictly positive real numbers and _M n (R) denotes the set 
of real square matrices of order n. 

Let us introduce the vector field F defined on ft = (Rl)™: 


F(x) 


/ /i(a;i,... ,x n ) 

\ fn(x 1 , . . • , Xn) 
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with: 


n n 

fi{x 1 ... ,x n ) = ati JJ x 9 y - Pi Xj‘1 , i = l...n 
3=1 i=i 

F is C 1 and therefore locally lipschitz on the open 17. Cauchy-Lipschitz theorem 
ensures the existence and unicity of a maximal solution of S(a, P,G, H) in 17, 
given any initial condition s(0) = x° G 17. 

This definition of S-systems with power-law differential equations is strongly 
linked with equations of chemical kinetics. As an example, if we consider the 
following chemical pathway: 

A + 2B^C^3D + E 

then the mass-action law applied to species C gives the equation: 

ft =kiab 2 ~k 2 d 3 e 

(capital letters design species and small letters design concentrations) 
Therefore in definition 2.1, coefficients a, and Pi are sometimes called kinetic 
rates while gtj and hij are called kinetic orders. 

S-systems are part of a broader formalism known as quasi-monomial (QM) 
systems (see [2]). An interesting result shows that QM systems can be expressed 
in the form of Lotka-Volterra quadratic systems (see [2] for details). 

2.2 Equilibrium points 

The study of the phase portrait of a S-system S(a,P,G,H) begins with the 
search for equilibrium points in 17. To find them, we have to solve the system: 

n n 

ai n x T = & n ? i = 1 • • • n (!) 

3 = 1 3 = 1 

In this paper we will use the following notation: 

Given a vector x £ (R^)" and a real square matrix A = (a*j) 4 - =1 n , we define 
the vector x A € (R!j_) n by: 


(x A )i = U , i = l...n 

3=1 

With this notation, we can express equation (1) as follows: 

x G ~ H = b 

where b is the vector {P\/a\,... ,Pn/ex n )- 

Taking the neperian logarithm, this equation leads to: 

(< G-H ) Ins = In b 
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(the logarithm is applied to all components of vector, i.e. In x is the n-dimensional 
vector (lnxi,..., In£.„)). Posing y = lnx, we are brought back to the resolution 
of a ?r-dimensional linear system in y. 

We have therefore the following proposition: 

Proposition 2.1 A S-system S(a, (3,G, H) has a unique equilibrium, x in Q 
(i.e. a positive equilibrium) if and only if the matrix (G — H) is invertible, x 
can be calculated by the formula: 

x = b^-V 1 (2) 

2.3 Stability analysis of the equilibrium 

The stability analysis of the equilibrium x uses the study of the spectrum of 
Jf(x) (the jacobian of F in x). The question we tackle here is to find some 
relationship between the stability of x and some properties of the matrix G — H. 
As a motivating example, let us consider the one-dimensional case. A one¬ 
dimensional S-system is expressed by a single differential equation: 

S(a,/3,g,h) : x = f(x) = ax 9 — (3x h 

where a,(3> 0 and g, h e R. The positive equilibrium x of ( S ) exists and is 
unique if and only if g — h ^ 0. In this case, an obvious calculation leads to: 

|^(z) = ax 9 ~ 1 (g - h) 

so the stability of x depends directly on the sign of g — h: it is asymptotically 
stable if g — h < 0 and unstable if g — h > 0, regardless of parameters a and (3. 

In the n-dimensional case, the stability depends on the sign of the real parts 
of the jacobian’s eigenvalues. Derivating the functions fi{x \..., x n ), we obtain, 
for i, j = 1... n: 

^ = (fT- n *k k ■ - h v ) ( 3 ) 

j j fc=i 

As in one-dimensional case, we thus obtain a formula that links the jacobian of 
F in x with the matrix G — H. However, it is not trivial to link the spectrum 
of Jf(x) with the spectrum of G — H. 

Let us recall here the definition of stability of matrices: 

Definition 2.2 A real square matrix A of order n is said to be stable (resp. 
semi-stable) if all its eigenvalues A i, i = 1 ...n, have a negative (resp. non 
positive) real part. 

We could hope that the stability of matrix G — H was sufficient to deduce the 
stability of x. However this is not true, as we can see in the following example. 
For n = 2, consider the S-system: 

/ C1 n . / x = 3 xy 2 - 2x a 

' \ y = 4 x 3 y 4 — x 5 y 3 
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we have: 


a = 


3 

4 




1 2 
3 4 


H = 


4 0 \ 

5 3 j 


The matrix G — H is equal to: 

G-H = 


-3 2 
-2 1 


Its characteristic polynomial is y(A) = (A + l) 2 so the matrix is stable. 

Since G — H is invertible, there is a unique equilibrium: x = (4p, 2|p).We can 
calculate the two eigenvalues Ai and X 2 of Jp{x). We find that Ai,A 2 > 0, 
implying that x is an unstable node. As a result, in spite of the stability of 
matrix G — H , the equilibrium x is unstable. 

The stability of G — H is therefore insufficient to deduce the stability of x. 
We need a stronger property known as sign stability (see [9],[11]). 

Definition 2.3 Two real square matrices of order n, A = (fljj)j - =1 n and 
B = {bij) i j =1 n , have the same sign pattern if: 

Vi, j = 1...n, sgn(aij) = sgn(bij) 


The function sgn is the classical signum function: 


Vx € R, sgn(x) 


+1 if x > 0 

0 if x = 0 

— 1 if x < 0 


Definition 2.4 A real square matrix A of order n is said to be sign stable (resp. 
sign semi-stable) if all the matrices that have the same sign pattern are stable 
(resp. semi-stable) in the sense of definition 2.2. 


In [9] we find a characterization of the sign semi-stability: 

Theorem 2.1 (Quirk-Ruppert-Maybee) 

A real square matrix A = s**? 71 semi-stable if and only if it satisfies 

the following three conditions: 

(%)Vi = l...n, an < 0 

(ii) Vi ^ j , aijaji < 0 

(in) for each sequence of k >3 distinct indices ii,... ,ik, 

We have. l)i(2) • • • ^i(k—l)i(k)^i(k)i(l) ^ 

(The third condition is equivalent to the fact that the directed graph associated 
to A admits no k-cycle for k > 3) 

With this notion, we can formulate the following proposition, which links the 
stability of the equilibrium x of a S-system with the sign semi-stability of matrix 
G-H: 
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Proposition 2.2 Let consider a n-dimensional S-system S{a, (3,G, H). We 
assume that G — H is invertible and we note x the unique positive equilibrium 
of (S). We also assume that x is hyperbolic (i. e. none of the eigenvalues of the 
jacobian of F inx have zero real part). 

If the matrix G — H is sign semi-stable (i.e. if it verifies the three conditions 
of theorem 2.1) then, regardless of parameters a and /3, the equilibrium x is 
asymptotically stable. 


PROOF. 

Let us note J the Jacobian of F in x and P the matrix G — H. The equation 3 
yields: 


dfi H 

^—(x) = —Pij 
OXj Xj 


with 7 i = a, rifc=i x 9 k ik ■ As 7 * > 0 and Xj > 0 for all i and j, matrices J and P 
have the same sign pattern. We can thus deduce that J is semi-stable and as x 
is supposed hyperbolic, it is asymptotically stable. 


□ 


Let us remark that the latter equation gives, in matricial notation: 

J = TPD- 1 

where T and D are diagonal matrices: 


r = 

( 71 


, D = 

( X\ 

.. ) 



In ) 


V 

X n J 


so sgn(det(J)) = sgn(det(P)). As we have supposed that P is invertible, we 
deduce that J is also invertible and does not have null eigenvalues. So we 
supposed the hyperbolicity of x in order to avoid imaginary eigenvalues of J. 
We can easily verify in the previous example that G — H is stable but not sign 
senri-stable (522 — /i 22 = 1 > 0 ). 

3 Local approximation of dynamical system us¬ 
ing S-systems 

In this part, we propose an algorithm for approaching the equilibria of a dy¬ 
namical system using S-systems. Simultaneously, we obtain a S-system that 
approximates the initial system around the equilibrium. 
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3.1 Monomial approximation of a positive vector field 

(see [17],[14],[16]). 

Let’s consider the positive vector field F : (R?j_) n —> (R^)”. 


F(x) 


/ fi(xi ,..., x n ) 
\ fni.%1 j • • • 5 %n) 


We will suppose F sufficiently smooth on (R(() n . 

Let us define the following change of variables: y = In a:, and express the loga¬ 
rithm of F(x) as a function G of the new variable y: 


lnF(x) =lnF(e y ) = G{y) 

The function G is sufficiently smooth on R™. Given any arbitrary point y° £ R™, 
let us write the Taylor expansion of g t (for i = 1... n) in the neighborhood of 
y° at the first order: 


V* = l...n, gi(y) = gi(y°) + - y^-^iy 0 ) + o (|| y - y° 

3 =1 Vj 


We introduce the functions g t (y) for i = 1... n: 


n r, 

V* = 1 • • • n, gi(y) = gi(y°) + ^2(yj - y°)-^-(y°) 

j-i 


and the functions /» = exp(gi(y)): 

fi{x) = f e*(») 


U=1 ayj 


= e 9i(v°) 




As y = In a; and _g,(y) = In/* (a;), we have: 


fi(x) 


n (5 

a=i 


eg; 


(v°) 


x l 
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and: 


% 

% 


(y) = 


d Vj 




1 


d 


fi{e v ) dy 3 
1 ...dfi 


fi(x) 


oVi . 


dxj 


(fi(e y )) 

(e y ) 


dfi 


fi(x) dxj 

Therefore, we have defined a vector field F = 


(x) 


F(x) 



(4) 


a fix 0 ) = fi(x°) 9ij 

with: { i =1 

n ( X °) = — 

9lj{X) ffiafi) dx 3 [x } 


(5) 


The basic idea is to use the monomial vector field F as an approximation of F 
in a neighborhood of x°. 


Definition 3.1 Let F be a smooth n-dimensional vector field, F : (R+)" —> 
(M+)" and x° any vector of (R^) n . We call S-approximation of F in x° the 
vector field F defined by equations (4) and (5). 


The following proposition is basic for what follows: 

Proposition 3.1 Let F be a positive vector field and F its S-approximation in 
x°. The following equalities hold: 

• F(x°) = F(x°) 


dfi 


dfi 


VM = l...n, ^(x°) = 

CJXj C/Xj 

(or, which is equivalent: Jf{x°) = Jp(x 0 )) 


3.2 Finding equilibria of a dynamical system 

We consider a n-dimensional dynamical system of the form: 

(S) x = V+(x) -V~(x) 

where x lies in (R+) n and V + , V~ are positive vector fields. V + , V~ : (R+) ra —> 
(R!j_) n . For i = 1... n, the term id (a;) is the production term of the variable 
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Xi and v~(x) the decay term of x, t . We propose an algorithm for finding an 
equilibrium point of (S) that lies in (R+) n . Meanwhile, we get a S-system that 
approximates the system (S) around this equilibrium. 

Given a point x° in (R^) n , we introduce the fields V + and V~ which are 
the S-approximations of the fields V + and V~ in x°. Let us consider the n- 
dimensional S-system: 

(S x o) x = V + {x) - V~(x) 


using (4) and (5), we obtain: 


where: 


and: 


(S x o) : ii = on JJ x 9 ? 3 - fa JJ Xj ' 3 , i = l...n 

3 =1 3 =1 





vf{x°) n^?) m 

3 = 1 


3=1 


dvt 


vj ( x°) d Xj 


„,0 

X 3 


8V: 


Hi (x°' 


dxj 


(x°) 

(x°) 


( 6 ) 


(7) 


If the matrix G — H is invertible, the system (S x o) admits a unique equilib¬ 
rium x eq G (R^) n : 

X e q = b^r 1 


with b = (fii/ai,fin/an). This point x eq depends on the initial point x° 
where we made our approximation. Let x 1 = x eq be the new initial point where 
we make our new S-approximation. The algorithm (1) computes the iteration 
of that process. 


3.3 Correctness of the algorithm 


Let’s describe the first iteration. 

Let £° € (R+) n . With formulae (6) and (7), we define the quantities ai(x°), 
Hi(x°), gij{x°) and hij(x°).They depend on the choice of the initial point x°. 
We assume that the constructed matrices G(x°) and H(x°) verify the condition: 
det(G — H) ^ 0. Thanks to this assumption, there exists a unique equilibrium 
point of the system (S x o). We will denote it x 1 , and we define the function 
’L : (R+)" —> (R^) n that, to each x° G (R^j_) n associates the point x 1 . 

Our algorithm is iterative, in the sense that it computes: 


(I) 


x° g (r;)” 

x n +i = 
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This iterative process converges towards fixed points of T. However we do not 
a priori know if all fixed points of T are indeed limits of (I). In other words, 
we must find which fixed points are attracting. 

The correctness of the algorithm (1) is a consequence of the two following lem¬ 
mas: 

Lemma 3.1 The equilibria of initial system ( S ) are the fixed points of the func¬ 
tion T 

Lemma 3.2 Given a fixed point x o/T, there exists some initial points x° that 
lead to x by the iteration (/). In other words, the positive equilibria of ( S ) are 
the attracting fixed points of d'. 

PROOF. 

(First lemma) Let x £ (Rl) n such that det(G(;r) — H{x)) is different from zero, 
(for convenience, we will omit the dependency in x , and note for instance G in 
place of G(x)). Using equation (2), we have: 


xS,(x) = b ( - G ~ H '> 1 


( 8 ) 


where b is the vector (/3j/ctj) i=1 n 
Therefore: 


^(x) = x 



GG-H )- 1 = x 
b = x^- H ^ 


Vi = 1... n, 
Vi = 1... n, 


ii 


— hij 


3 =1 


»iK“=«-iur 

j =i t=i 


By definition, o ,; x 9 G (resp. /3* ]~[y=i x> j >3 ) the S-approximation of V + 
(resp. V~) in x. Proposition 3.1 implies then: 


T (x) = x 


V + (x) = V {x) 


Thus, the equilibria of (5) are the fixed points of the function U/. 


□ 

In order to prove the second lemma, we will use the following fixed point 
criterion: 

If the function ’F is a contraction on the open set W and if x £ W is a fixed 
point of H/, then x is the unique fixed point of T in IF and it is attracting, that 
is to say, for all x° £ W, the iteration (J) converges towards x. PROOF. 
(Second lemma) Let a; be a fixed point of T. We assume that det (G(x)-H(x)) 

0. The continuity of the determinant implies that there exists a neighboorhood 
W of x in which det(G — H) 0. To prove that x is attracting, it is sufficient 
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to show that 'L is contracting in a neighboorhood of x. For that, we show that 
the jacobian of 'F in x is zero. 

Using (6) and (7) and posing: 


' U+ = log(V+) 

U~ = log(U-) 

u = U+-U-= log 

P = G-H 


we obtain, for all x G W: 


'Pi(x) 


n 


i=i 



„(-!) 


(x) 


where [p\- 1 M is the inverse of the matrix P = G 

Let’s calculate p,-, 11 : 

* L J 


= G - H. 


Pij 9ij kij 

1 dvf 1 dvj 


V■ dxj V. 

(duj du~ 

Xj V dxj d Xj 

dui 

= Xj -— 

dxj 

in matricial notation: P = J u (x)A where J u (x) is the jacobian of the function 
U evaluated in x and A is the diagonal matrix: 




A = 


/ 


V 


Xi 


Therefore P 1 = A 1 ( J u (x )) 1 = A 1 (J u -i(x)) ( u 1 is the reciprocal func¬ 
tion of u) and so: 

w- ■ i (-i) 1 du ~ l 

\h, 3 = 1 ...n, Pij = 


Xi dxj 


with (8) we have, for i = 1... n and x £ W: 


&i(x) = Xi exp ( - — Uj ( x ) (x) 


j=i 


Xi " ' ' dxj 

,-i 




i=i 


Xi dxi 
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Deriving this (and omitting the dependency in x), we get, for k ^ i: 

Uj(x) du - 1 


dxk 


E 1 

3=1 


d 2 u „• 1 


dxjdxk 


ex p - E 


3=1 


dxj 


and 


d*j 

dxi 


E 


d 2 u • 


du- 1 

^ '~ J dxjdxk Xi Uj dx.j 
1 — 1 J 1 = 1 J 


exp 


-E 


3=1 


Uj{x) du 
Xi dx-j 


-l 


As we have shown that the fixed points of \I/ are the equilibria of (5), we deduce 
that Vfc = 1... n, Uk{x) = 0, therefore: 

Jy(x) = 0 


We deduce that 'h is contracting in a neighboorhood of x, and then that x is 
attracting. This concludes the proof of the second lemma and the correctness 
of the algorithm. 

□ 


3.4 An example with multiple positive equilibria 


We present here the application of our algorithm for a dynamical system having 
multiple positive equilibrium points. It is a system known as biological switch 
(see [3]). 

Let’s consider the two dimensional dynamical system: 


x = 


V = 


— x 


1 + y 2 

6.75 

3.375 + a; 3 


(9) 


It represents the temporal evolution of two positive quantities x and y with 
linear decay and sigmoidal production (we use here the Hill function H~(z) = 
K n 

— - often used by biologists to model sigmoidal interactions). As we can 

K n + z n 

see on figure 1, This system shows three equilibrium points. The values of these 
points can be calculated: 


PI 


0.697 \ 

1.818 J 



P 3 


2.802 \ 
0.266 J 


We can show that P2 is unstable whereas PI and P3 are stable (cf. [3]). 

Applying our program in Maple, we found three different initial conditions 
each of which tends towards one of the three equilibrium points (see figure 1 
and numerical results below). The convergence appears to be fast since we need 
only 4 iterations to approach the equilibria with a precision of 10” 5 . We will 
discuss about the convergence speed in part 5.1. 
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X 


Figure 1: nullclines of system (9). (The dash line represents fi(x,y) = 0 and 
the solid one represents f'iix, y) =0). The central equilibrium P2 can be shown 
to be unstable while the two others, PI and P 3 are stable. The arrows represent 
the three experimentations described in the text. 


• With initial condition a; 0 = (2,2), algorithm finished in 4 iterations and 
found P2 with a precision of 10 -5 . The numerical S-system obtained is 
given by: 

( x = 1.500y _1 — x 
1 y = 1.837 x- 1 - 5 -y 




With initial condition x° = (0.2,1.5), algorithm finished in 4 iterations 
and found PI with a precision of 10 -5 . The numerical S-system obtained 
is given by: 

x = 1.745 y -1 ' 535 — x 
y = 1.647 x~ 0 - 274 — y 




With initial condition x° = (2, 0.2), algorithm finished in 4 iterations and 
found P3 with a precision of 10~ 5 . The numerical S-system obtained is 
given by: 

f x = 2.352 y-° 132 ~x 
\y = 3.879 x~ 2G ~y 
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3.5 Stability analysis of approximate S-systems 


Consider the n-dimensional dynamical system: 

( ^ = F(x) = V + (x)-V~(x) (10) 

1 x g (k;)" 

Algorithm 1 ensures that, given any initial condition x° in (R(j_) n , unless we fall 
in a degenerate case, we produce a sequence (x q ) q& fq (with x q G (R+)") that 
tends towards a limit point x G (Rl)" which is an equilibrium of (10). More 
precisely, x q = ^/^^(a; 0 ). 

Meanwhile, at each step, it provides us with a S-system S q {a q , P q , G q , H q ) which 
comes from the S-approximations of functions V + and V~ in x q . Thus, we have: 


OLq 

= a(x q ) 


Pq 

= P(x q ) 


Gq 


with: 4 = 

Hq 

II 

Co. ' 

II 

with: /if. = 
..n l J 


9ij(x q ) 

M x q ) 


where a, p, gij and hij are the functions defined in (5). If we assume that V + 
and V~ are at least C 1 . we deduce that these sequences converge, as q tends to 
oo, towards: 



—> a(x) 

def 

_ 

CXq 

— 

a 

P q 

-> P(x) 

def 

p 

G, 

- G(x) 

def 

G 

. H, 

-> H (x) 

def 

H 


Let ( S ) be the following S-system: 

n n 

( S ) : Xi = on xf° - % ]^[ x , * = 1 ... n ( 11 ) 

i=i j =i 


We want to know in what sense the system (11) approach the system (10). An 
answer is given by the following proposition: 

Proposition 3.2 F is supposed C r (r > 1). The equilibrium x of (10) is an 
equilibrium of (11). Moreover, ifx is hyperbolic, then the flow generated by (11) 
is topologically conjugate to the flow generated by (10) in a neighborhood ofx. 


PROOF. 

The first assertion is obvious with proposition 3.1. Let prove the second as¬ 
sertion: it is a direct consequence of the Hartman-Grobman theorem (see for 
instance [18]). 

Proposition 3.1 shows that systems (10) and (11) have the same linearized dy¬ 
namical systems in x. Thanks to the Hartman-Grobman theorem, we know that 
these systems are topologically conjugate to their linearized dynamical systems. 
By transitivity of the topological conjugation, we deduce that (10) and (11) are 
topologically conjugate around x. 
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□ 


This proposition implies that the stability of x for system (11) is the same that 
the stability of x for system (10). As an exemple, let us consider the following 
2 -dimensional dynamical system: 

2 4 

x x y 

2 + y (3 + x)(4 + y 3 ) 

5x 2xy 3 

3lx (x T l)(y + 2) 

We find the equilibrium point x « (1.2301,1.6950) and the matrix G — H: 

-0.709 -2.812 \ 

0.261 -2.541 ) 

Thanks to theorem 2.1, we see that G — H is sign semi-stable. The point x, as 
equilibrium of (Ex) is hence stable. 



4 Application to a biological example 

We present here a current work we are doing in collaboration with G. Curien (see 
[4]). The goal of this work is to understand the metabolic system responsible for 
the synthesis of aminoacids in Arabidopsis Thaliana. So far, we have focused 
our study on a subsystem of 14 variables, with 9 symbolic parameters. The 
differential equations present several strongly nonlinear terms due to allosteric 
control of some enzymes ; in particular, Hill functions and compositions of Hill 
functions. Since the latter are rational functions, seeking positive equilibria is 
equivalent to solving a polynomial system. Algebraic manipulations have led 
us to a simplified system with 5 polynomial equations in 5 variables. Because 
of the complexity of these equations, we were not able to achieve the resolution 
of this system with purely symbolic computations and manipulating symbolic 
parameters (we used Maple 8 and 9). That is why symbolic-numeric methods 
appeared as a satisfactory way to tackle this problem. As it is a system of equa¬ 
tions coming from biochemical kinetics, S-systems seemed to be an appropriate 
tool in this work. 

In vivo , this system exhibits a stationnary behavior. Giving realistic values of 
parameters, we managed, thanks to our algorithm, to find this positive equi¬ 
librium. We now have to study the S-approximation of the system near this 
equilibrium, with different realistic sets of parameters. An interesting idea is 
also to propose a piecewise S-approximation of the system in order to reproduce 
its behavior in a wider zone of the phase space. This work is in progress. 
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5 Discussions and concluding remarks 

5.1 Convergence of our algorithm 

The algorithm described above computes the iterations of a vectorial function 
\I/ on an initial point x° £ (R+) n , in order to converge towards a fixed point 
of \I/. As the jacobian of HI is the null matrix in those fixed points, we know 
that the convergence speed is very fast (up to four or five iterations in all the 
examples presented, for a precision of 10 -4 or 10~ 5 ). As a matter of fact, we 
are in a case where the speed of convergence is the best possible. Indeed, if the 
function <3/ is AT-contractant, one can easily verify that the convergence of the 
iteration is in K n (where n is the number of iterations). Since J^(:r) = 0, then 
we can find a neighborhood of x wherein T is A'-contractant for any 0 < K < 1. 
However, even if the speed of convergence is very fast, the algorithm behaviour 
is strongly dependent on the choice of initial point x°. Indeed, if initial system 
has multiple positive equilibria, each of them have distinct basins of attraction. 
We cannot a priori know in which of these basins is the point x°. We even 
cannot ensure that x° actually lye in one of them. In fact, the study of basins of 
attractions of such iterations is a complex issue. The boundaries of such basins 
can be quite complicated, even fractals [8]. As an example, we launched our 
algorithm for the switch system (equations (9)) with initial conditions taken on 
a grid of ]0,4] 2 . To vizualize the three basins, we colored the initial points (fig 
2 ). 

5.2 interaction between symbolic and numerical calculus 

As we said in the introduction, a large part of research concerning the analy¬ 
sis of biological phenomena uses both symbolic and numerical techniques. The 
S-systems as we described represent a large class of systems, yet their simple 
mathematical expression allows symbolic manipulations, providing a practical 
framework of study. Algorithm 1, as presented here needs numerical estimations 
of symbolic parameters. Nevertheless the technique of S-approximation (def 3.1) 
consists of symbolic manipulations (in particular, we use symbolic computation 
of partial derivatives). It can be calculated in any point of the phase space and 
can include symbolic parameters. 

S-approximation gives a computable and rather good approximation of ODE 
systems (see [17] for a comparison between power-law approximation and lin¬ 
earization). A very interesting idea is therefore to use the context information 
(given for instance by biologists) of a particular system in order to create a 
piecewise S-approximation of this system. This should provide a global approx¬ 
imation interpolating the system in some critical points in the phase space (see 
[13]). 


16 




Figure 2: Basins of attraction of points PI (dark), P 2 (white) and P3 (grey). 
We obtained these graphs by applying algorithm 1 for system (9) with initial 
conditions taken in a regular grid of ]0,4]. 
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Algorithm 1 Search of an equilibrium point of system ( S ) 

Require: 


X = x° £ (R;!j_) n (initial condition) 

V + , V~ : positive vector fields defined over (R^_) n 
e > 0 : precision 

Ensure: unless we fall in a degenerate case, we find a point y close to a positive 
equilibrium of ( S ) with the precision e. Meanwhile, we obtain the S-system 
( S y ) that approximate system ( S ) around this equilibrium. 

repeat 

Y := X 

for i = 1 to n do 
for j = 1 to n do 

Xj dv+ 

9ij : v+(X) dX 3 

h — Xj dVi 
v~{X)dX J 

end for 

n 

a i ■■= v+{X) l[(Xi)-^ 
j=1 

n 

Pi ■■= v~(X) U(X^ 
i=i 

bi '■= Pi/a-i 

end for 

if det(G — H) ^ 0 then 
X := b^r 1 

else 

degenerate case: algorithm terminated —► restart the algorithm with a 
new initial condition 

end if 

until || X — Y ||< e 
Result := X 





